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Summary 


The  evolution  of  unsteady  boundary  layers  in  the  vicinity  of  the  leading  edge 
of  a  thin  oscillating  airfoil  has  been  examined  with  a  novel  numerical  method 
which  is  able  to  deal  with  the  movement  of  the  stagnation  point  and  with  reg¬ 
ions  of  reverse  and  separated  flow,  solutions  to  the  unsteady  boundary- layer 
equations,  with  a  prescribed  pressure  distribution  which  causes  flow  reversal 
and  separation,  demonstrate  the  importance  of  numerical  steps  in  distance  and 
time  and  that  a  requirement  similar  to  the  stability  criterion  of  Courant, 
Friedrichs  and  Lewy  must  be  satisfied  to  avoid  numerical  errors.  At  the  lower 
reduced  frequencies  of  the  investigation,  solutions  could  not  be  obtained  with 
this  procedure  and  it  was  necessary  to  introduce  Interaction  between  the  vis¬ 
cous  and  lnvlscid  flows.  The  solutions  obtained  with  the  interactive  method 
were  Increasingly  different  from  those  without  interaction  as  the  reduced  fre¬ 
quency  was  decreased  towards  zero  and,  for  some  combinations  of  Reynolds  number 
and  frequency,  exhibited  behavior  consistent  with  the  instability  of  separation 
bubbles. 
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1.0  INTRODUCTION 


The  lift  and  drag  characteristics  of  airfoils  at  moderate  Reynolds  numbers  can 
be  affected  by  separation  bubbles  which  occur  close  to  the  leading  edge  and, 
at  high  angles  of  attack,  can  increase  in  size  to  cause  stall.  The  added  com¬ 
plexity  of  unsteady  motion  such  as  that  associated  with  the  rotor  blades  of 
helicopters  implies  that  the  flow  characteristics  are  influenced  by  amplitude 
and  frequency  and  that,  in  particular,  the  stall  characteristics  can  be  consid¬ 
erably  modified.  The  investigations  of  Carr,  McAlister  and  McCroskey  (1977), 
Francis,  Keese  and  Retelle  (1983),  Daley  and  Jumper  (1984)  and  Lorber  and 
Covert  (1986)  examined  these  effects  over  limited  ranges  of  the  parameters  and 
that  of  Carr  et  al.  (1977)  provides  detailed  information  of  the  mechanism  of 
dynamic  stall  of  an  oscillating  airfoil.  It  appears  that  stall  is  associated 
with  flow  reversals  in  the  unsteady  boundary  layer  and  that  these  may  translate 
downstream  or  upstream  depending  upon  various  parameters  including  the  radius 
of  the  leading  edge  of  the  airfoil.  At  some  stage  in  the  cycle,  stall  occurs 
and  is  preceded  by  a  vortex  which  forms  close  to  the  surface  and  is  probably 
associated  with  a  breakdown  of  the  unsteady  boundary  layer. 

The  above  physical  problems  involve  laminar,  transitional  and  turbulent  flow 
and  their  representation  requires  a  numerical  calculation  procedure  which  can 
provide  accurate  solutions  to  conservation  equations  in  all  regions  of  flow  as 
well  appropriate  transition  and  turbulence  models.  Here  we  are  concerned  with 
the  numerical  solution  procedure,  its  development  to  represent  the  regions  of 
reverse  flow  and  use  to  examine  the  nature  of  solutions  for  parameters  close 
to  those  associated  with  stall.  The  emphasis  is  on  regions  of  flow  close  to 
the  leading  edge  of  a  thin  oscillating  airfoil  and  calculations  are  performed 
with  prescribed  pressure  gradient  and  with  interaction  between  solutions  of  the 
viscous  and  lnvlscid  equations.  With  the  configuration  chosen,  an  analytical 
solution  for  the  potential  flow  equations  was  available. 

Previous  consideration  of  steady  boundary  layers  and  their  solution  by  an 
interactive  procedure,  has  been  reported  by  Cebeci,  Stewartson  and  williams 
(1981)  for  a  model  problem  consisting  of  a  thin  ellipse  at  incidence.  Their 
study  showed  that  the  solutions  were  well  behaved  and  unseparated  provided  the 
angle  of  attack  was  less  than  1.155  degrees.  At  higher  angles,  separation 
occurred  with  an  associated  singularity  which  was  overcome  by  the  use  of  the 


Interactive  procedure  and  results  were  obtained  Cor  small  regions  of  separated 
Clow.  There  is,  however,  a  limiting  size  oC  separation  bubble  beyond  which 
Cebeci,  et  al.  (1981)  could  not  obtain  solutions  and  this  may  be  related  to  the 
physical  phenomenon  oC  open  separation  and  stall.  A  similar  result  was 
obtained  by  Stewartson,  Smith  and  Kaups  (1982)  who  used  a  triple-deck  approach 
and  Cound  that  their  calculations  oC  separation  bubbles  could  break  down  with 
a  small  increment  in  pressure  gradient.  They  also  observed  that  their  solu¬ 
tions  were  not  unique  and  their  results  may  imply  that  large  separation  bubbles 
cannot  exist  in  laminar  Clows  at  high  Reynolds  numbers. 

The  unsteady-Clow  calculations  reported  here  were  obtained  with  Keller's  box 
method  (1974)  Cor  the  solution  oC  the  boundary- layer  equations.  In  regions  oC 
Clow  reversal,  a  requirement  similar  to  the  stability  criterion  oC  Courant, 
Friedrichs  and  Lewy  (CPL),  see  Isaacson  and  Keller  (1966),  is  satisCied  by  the 
use  oC  the  characteristic  box  procedure  discussed  by  Keller  (1978)  and  Cebeci 
(1986)  and  the  interactive  procedure  is  based  on  the  Hilbert-integral  previ¬ 
ously  used  by  Cebeci  et  al.  (1981).  This  combination  oC  methods  represents  the 
best  possible  approach  available  to  the  authors  and  allows  the  importance  oC 
the  stability  criterion  to  be  examined  as  well  as  the  structure  oC  the  solu¬ 
tions.  OC  necessity,  a  limited  number  oC  parameters  is  considered  but  encom¬ 
passes  a  range  oC  relevance  to  oscillating  airCoils. 

The  Clow  conCiguration  under  consideration,  equations,  initial  and  boundary 
conditions  are  examined  in  the  Col lowing  section  which  is  Collowed  by  a  brieC 
description  oC  the  solution  procedure.  The  results  are  presented  and  discussed 
in  Section  4  and  the  paper  ends  with  a  summary  oC  the  more  important 
conclusions. 
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2.0  FLOW  CONFIGURATION,  CONSERVATION  EQUATIONS, 
INITIAL  AND  BOUNDARY  CONDITIONS 


We  consider  flow  over  an  ellipse  with  a  thickness  ratio  x(=b/a)  much  less  than 
unity  at  an  angle  of  attack  a.  The  body  surface  is  defined  by 

x  =  -a  cos0,  y  =  ax  sin0  -*  <  0  <  t 

and  the  corresponding  external  velocity  for  steady  flow  can  be  deduced  from 
inviscid  flow  theory  to  be 


u°(s,t) 

e 


ilia. 

/i  +  s2 


(i) 


o  o 

Here  ue(s,t)  denotes  a  dimensionless  velocity,  u  /u^d  +  x),  the  parameter 

£  denotes  a  dimensionless  distance  from  the  nose  related  to  the  x-  and 

2  2 

y-coordinates  of  the  ellipse  by  x  =  1/2  ax  £,  y  =  ax  £,  and  (=a/x) 
represents  a  dimensionless  angle  of  attack.  The  parameter  l  is  related  to 
the  surface  distance  s  by 


s 


ax 


I1  (1  +  l2)1/2dl 
0 


The  boundary- layer  equations  for  unsteady  incompressible  laminar  flows  on 
oscillating  airfoils  are  well  known  and  can  be  written  as 


3u  3v 
3s  +  3n 


(2) 


3u  A  3u  ^  3u 
r  +  u  r  +  v  r 
3t  3s  3n 


3u  3u  ,2 

_ e  _ e  3  u 

3t  ue  3s  V  ,  2 
3n 


(3) 


Solutions  to  these  equations  are  usually  obtained  for  prescribed  boundary 
conditions  given  by 


n*0,  u  =  v  *  0;  n=n,  u*u  (s,t)  (4) 

e  e 

and  we  shall  refer  to  this  as  the  standard  problem.  In  the  interactive  problem 
we  determine  u  (s,t)  partly  from  inviscid  theory  and  partly  from  the  pressure 


distribution  resulting  from  the  blowing  velocity  d/ds  induced  by  the 

boundary  layer.  Thus  we  write 


u  (s,t)  =  u  (s,t)  +  u  (s,t) 
e  e  c 

where  u°(s,t)  is  the  inviscid  velocity  and  u  (s,t)  is  related  to  the 
©  c 

blowing  velocity  by  a  variation  of  the  Hilbert  integral 


u(s,t)  »  J  I  (u  i‘) 
c  »  _  as  e  s-ff 


(5) 


(6) 


which  is  valid  for  straight  walls  but  can  be  generalized  to  airfoils  as  dis¬ 
cussed  by  Cebeci  and  Clark  (1984).  The  freestream  velocity,  consistent  with 
Bq.  (1)  has  the  form 


u°(s,t) 

© 


5+5  (1  +  A  sinwt  ) 

_ o _ 


/l  +  52 

where  A  is  an  amplitude  parameter  and  u  is  a  dimensional  frequency. 


(7) 


For  attached  flow,  the  effect  of  uc(s,t)  is  generally  weak  but  is  enhanced 
in  the  neighborhood  of  separation  as  can  be  surmised  by  noting  that  the  inte¬ 
grand  in  Bq.  (6)  would  otherwise  develop  a  strong  singularity  at  separation  and 
cause  the  solutions  to  break  down  further  downstream.  As  discussed  by  Cebeci, 
et  al.  (1981),  it  is  sufficient  to  replace  Bq.  (6)  by 

s  (u  «*)' 

uc(s't)  -  i  /  r^~r dB  (8) 

sa 

where  the  prime  denotes  differentiation  with  respect  to  s  and  s  and  s. 

a  d 

denote  the  beginning  and  the  end  of  the  interaction  region. 


To  complete  the  formulation  of  the  problem,  upstream  boundary  conditions  must 
be  specified  in  the  (t.n)  plane  at  some  s  ■  sq  as  well  as  Initial  conditions 
in  the  plane  (s,n)  at  t  a  0.  If  steady-flow  conditions  prevail  at  t  =  0,  the 
initial  conditions  can  be  obtained  easily  for  both  surfaces  by  solving  the 
conservation  equations  for  steady  flow  which,  in  this  case,  are  given  by  Eq. 
(2)  and  by 
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There  is  no  problem  with  the  initial  conditions  for  Eqs.  (2)  and  (9)  since  the 
calculations  start  at  the  stagnation  point  where  they  admit  similarity 
solutions. 

The  generation  of  the  upstream  boundary  conditions  for  Eqs.  (5)  and  (6) 
requires  a  special  numerical  procedure,  since  the  complete  velocity  profile 
distribution  on  a  previous  time  line  is  known,  solutions  can  be  determined  on 
the  next  time  line  by  an  explicit  method.  If  we  wish  to  avoid  stability  prob¬ 
lems,  however,  an  implicit  method  is  required  and  generation  of  a  starting 
profile  on  the  new  time  line  becomes  a  problem. 

In  order  to  explain  the  problem  further,  it  is  instructive  to  see  what  happens 
to  the  stagnation  point  as  a  function  of  time.  For  this  purpose  let  us  con- 

o 

slder  Eq.  (7).  Since  ufi  =  0  at  the  stagnation  point,  its  location,  £s»  based 
on  the  external  streamlines  is  given  by 

l  -  -%  (1  +  A  sinwt)  (10) 

s  o 

Figure  1  shows  the  variation  of  the  stagnation  point  with  time  for  one  cycle 

according  to  Eq.  (10)  with  A  =  1,  u  =  x/4.  We  see  that  when  t  =  2,  the 

stagnation  point  £  is  at  -2\  ,  and  when  t  =  6,  it  is  at  0,  etc.  If  £ 
SO  s 

were  fixed,  we  could  assume  that  u  =  0  at  5  =  for  all  time  and  all  n,  but 

o 

this  is  not  the  case.  It  is  also  possible  to  assume  that  the  stagnation  point 
is  coincident  with  zero  u-velocity  for  a  prescribed  time  but  we  should  note 
that  the  stagnation  point  defined  by  Eq.  (10)  is  based  on  the  vanishing  of  the 
external  velocity.  For  a  time-dependent  flow,  this  does  not  imply  that  the 
u-velocity  must  be  zero  across  the  layer  at  a  given  ^-location  and  specified 
time;  Indeed  flow  reversals  can  occur  due  to  the  movement  of  the  locus  of  zero 
u-velocity  across  the  layer  and  can  cause  numerical  instabilities  which  require 
the  use  of  special  numerical  schemes  as  discussed  by  Cebeci  and  Carr  (1981). 


3 . 0  SOLUTION  PROCEDURE 

With  the  upstream  boundary  conditions  determined  by  the  procedure  of  Cebeci  and 
Carr  (1981)  and  with  the  initial  conditions  obtained  from  the  solution  of  Eqs. 
(2)  and  (9)  subject  to  the  boundary  conditions  given  by  Eqs.  (4)  to  (6),  Eqs. 
(2)  and  (3)  can  be  solved  for  both  standard  or  inverse  problems.  In  practice 
a  standard  procedure  is  used  up  to  a  specified  ^“location  after  which  the 
calculations  may  proceed  by  either  standard  or  inverse  procedures.  For  exam¬ 
ple,  the  evolution  of  the  boundary  layer  on  an  oscillating  airfoil  with  pre¬ 
scribed  pressure  distribution  is  determined  with  the  standard  procedure  and  the 
Inverse  procedure  is  used  after  a  short  distance  from  the  leading-edge  region 
where  the  inviscid  and  viscous  flow  equations  are  solved  interactively. 


To  solve  the  equations  for  both  standard  and  inverse  problems  we  use  modified 
forms  of  Keller's  box  scheme.  The  Mechul  function  formulation  of  Cebeci  (1976) 
is  used  in  the  inverse  case  and  treats  the  external  velocity  as  an  unknown. 
Before  we  describe  this  formulation  and  the  solution  procedure,  it  is  conveni¬ 
ent  to  write  Eqs.  (2)  to  (4)  in  a  form  more  suitable  for  computation.  To 
achieve  this  we  define  dimensionless  distances  n  and  I  and  time  x  by. 


n  *  [ 


R( 1  +  t.)  1/2 
- L_]  Q 

2t? 


u  (1  +  t  ) 

o>  1 


(11a) 


with  R  =  2au  /v,  and  a  dimensionless  stream  function  f  by 

OB 

2  1/2 

♦(s,n,t)  =  ((1  +  t.)au  vt.]  f(s,n,T)  (lib) 

i  “l 

These  relations  may  be  introduced  into  Eqs.  (2)  and  (3)  to  give,  with  primes 
denote  differentiation  with  respect  to  n, 


where 


A  3 w  ^  3w 
f " '  +  r-  +  w  —  = 
3t 

3s 


a  f ' 

3x 


+  f' 


3s  3s 


w 


u 

e 


V1  +  v  • 


f  • 


“J1  +  v 


(12) 


The  boundary  conditions  follow  from  Eqs.  (4)  to  (6)  and  can  be  written  as 


n  -  0,  f  =  f’  =  0, 


(13) 
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M 

•m 

M 


$$9 

wig 


n  =  n  ,  f'  =  w, 
e 


-o  .  ,  b  dA  da 

w  =  u  +  e  J - 

Q  —  —  —  — 

s  ds  s  -  a 
a 


Here  A  denotes  a  dimensionless  displacement  thickness  given  by 

A(s,t)  =  n  w  -  f  (1 

e  e 

and  e  is  a  parameter  defined  by 

1  r  2  .1/2 

C  ~  tt1  R( 1  +  tL)]  U 

To  pave  the  way  for  the  description  of  the  numerical  method,  we  define  a  new 
variable 

0=  *  a 


and  write  Eq.  (12)  as 

3w  3w  3f '  3f ' 

f"'  +  f'0  +  r  +  w  r  =  t —  +  f’  t —  (17) 

3t  3s  3t  3s 

where  the  overbar  has  been  omitted.  We  use  Keller’s  box  method  (1974)  to  solve 
this  equation.  In  regions  of  no  flow  reversal  the  so-called  standard  box 
method  is  used  and  where  there  is  flow  reversal,  this  is  replaced  by  the  char¬ 
acteristic  scheme  which  is  based  on  the  solution  of  Eq.  (17)  along  streamlines 
as  described  by  Keller  (1978)  and  Cebeci  (1986).  This  scheme  allows  the  step 
sizes  in  the  t  and  s-directions  to  be  automatically  adjusted  to  ensure  that 
the  region  of  backflow  determined  by  the  local  streamlines  does  not  violate  a 
condition  like  the  Courant,  Friedrichs  and  Lewy  (CFL)  stability  criterion. 
Although  the  zig-zag  scheme  of  Krause  et  al.  (1968)  can  also  be  used  for  this 
purpose,  it  can  be  inaccurate  in  regions  of  large  flow  reversal  since  the  ori¬ 
entation  of  numerical  mesh  is  chosen  a  priori,  as  discussed  by  Cebeci  (1986) 
and  later  in  this  paper. 

To  solve  Eqs.  (17)  and  (14)  with  the  box  scheme  and  the  Mechul  function  formu¬ 
lation,  we  let 


f'  =  e 


and  Introduce  a  new  function  g  defined  by 


e '  =  g 


(19a) 


and  with  w(x)  treated  as  unknown 


w*  =  0 


(19b) 


and  write  Eqs.  (16)  and  (17)  and  their  boundary  conditions  as 

3e 


0'  = 

g*  +  g0  + 
n  =  0,  f  =  e  =  0; 


3s 

3w  3w  3e 


3e 


3t  +  w  3s  3t  +  e  3s 


(19c) 

( 19d) 

(20) 


..  o  .  ,  dA  da 

"  ’  V  «  -  w.  “  -  ‘  J  £ 

sa 

To  write  the  difference  equations  for  the  system  given  by  Eqs.  (19)  and  (20), 
we  consider  a  net  cube  in  which  the  net  points  are  denoted  by 


-  0, 

s, 

3  s,  , 

+ 

r 

i  - 

1. 

2,  .  .  .  r  I 

i 

i-1 

i 

=  o. 

T 

*  T 

+ 

k 

n  * 

1. 

2  f  •  •  •  t  N 

n 

n-1 

n 

*  0, 

"  Vi 

+ 

hj 

j  - 

1. 

2 1  •  •  •  f  J 

=*  As. , 

k  *  At 

and  h . 

3 

An .. 

i 

n  n 

j 

j 

(21) 


The  difference  approximations  to  represent  Eqs.  (19a)  and  (19b)  are  obtained  by 

averaging  about  the  midpoint  (s.,  t  ,  n  ), 

i  n  j~l/2 


-1.  l,n  i,n. 
hj  ej  '  ej-l} 


i.n 

9j-l/2 


.  "I#  i#n  if n. 
h3  (Wj  -  wj_1)  -  0 


(22a) 

(22b) 
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where,  for  example, 


i,n  1  ,  i,n  i.n 
ej-l/2  ’  2(‘}  *  'l-l> 


The  finite  difference  approximations  to  Eqs.  (19c)  and  (19d)  are  obtained  by 

centering  all  quantities  except  9  at  the  center  of  the  cube  (s^_  ^  ,  Tn-l/2' 
n^  )  by  taking  the  values  of  each  say  q,  at  the  four  corners  of  the  box,  that  is. 


i-l/2,n  ,  i.n  i-l,n.  ,  i.n  i-l,n. 

qj-l/2  =  \  (qj-l/2  +  q J-1/2*  (qj  +  qj-l  +  qj-l  * 


(24a) 


The  centering  of  0  is  done  by  writing  it  as 


Qi-l/2,n-l/2  ,  ,_i-l/2,n  .i-l/2,n-l/2x 

3-1/2  -  ^  ‘  3  *  >1  1 


(24b) 


In  this  notation,  the  difference  approximations  to  Eqs.  (19c)  and  (19d)  can  be 
written  in  the  form: 


hl  (9j  9j-i>  =  Ci  (6i  ei*l) 


(25a) 


hj  (gj  ~  9j-l)  +  gJ-l/20j-l/2  +  kn  (wn  wn-l)  +  ri  twj-l/2(wl  wi-l)] 

3  kn1(;n  '  Vl>  +  rilt«j-l/2<«i  -  Vl>  "  qj-l/2(?i  "  ?i-l)]  (25b) 

where,  for  example. 

-  a  i-l/2,n-l/2  ~  _  i-l/2,n  “  _  i,n-l/2  A  _  fli-l/2,n-l/2  .... 

ej  *j  '  en  -  e j-1/2  '  ei  ~  j-1/2  '  0j  -  8j  (26) 

Following  the  procedure  of  Cebeci,  et  al.  (1981)  the  boundary  condition  involv¬ 
ing  the  Hilbert  Integral  in  Eq.  (20)  can  be  written  in  the  form 

i.n  .  i,n  _i,n.  i,n  . _  . 

WJ  *  cii  e(nJ  WJ  ’  CJ  }  =  TJ  (27) 

where  e  is  the  matrix  of  interaction  coefficients  defining  the  relationship 

between  the  dimensionless  displacement  thickness  and  external  flow  and  the 
i.n 

parameter  T_  represents  terras  where  values  are  assumed  to  be  known.  It 

J 

is  given  by 


—  i-1  1 

„i,n  ,  Ovi.n  .  r*  _  . m , n  .  .  m ,  n 

T2  -«V  *  ‘  Cl"  ‘Jui  ^ 


6343H 


9 


To  compute  the  additional  unknown  of  Eq.  (27),  we  write  Eq.  (10)  in  the  form 


i.n 

*3-1/2 


(29) 


so  that  the  system  consisting  of  Bqs.  (22),  (25),  and  (29)  can  be  solved 
subject  to  the  boundary  conditions  given  by  Eq.  (27)  and  which  follow  from  Eq. 
(20). 


f  =0  *  e  -  0; 

o  o  o 


e  a  w 


J  J 


(30) 


The  above  system  can  be  linearized  by  Newton's  method  and  the  resulting  linear 
system  solved  by  the  block  elimination  procedure  described  in  Cebeci  and 
Bradshaw  (1984)  for  both  standard  and  Inverse  formulations.  In  the  former 

o 

case,  it  is  sufficient  to  set  c  *  0  in  Eq.  (27)  so  that  w  is  equal  to  uft. 


Ve  follow  the  above  solution  procedure  when  there  is  no  flow  reversal  across 
the  layer.  If  separation  is  identified  from  the  values  of  u*'n,  use 
the  characteristic  scheme  which  has  recently  been  described  by  Cebeci  (1986) 
in  relation  to  the  standard  problem  of  computing  the  impulsively  started  lam¬ 
inar  flow  over  a  circular  cylinder.  The  solution  procedure  in  this  case  is 
similar  with  small  adjustments  due  to  the  manner  in  which  the  difference  equa¬ 
tions  are  adjusted  to  the  modified  form  of  Eq.  (19d).  Noting  the  definition 
of  local  streamlines,  we  write 


,  ds 

dT  3  - 

e 


(31) 


If  we  denote  distance  in  this  direction  by  q  and  the  angle  that  it  makes  with 
the  T-axls  by  a,  then  Eq.  (19d)  can  be  written  as 


g'+g8  +  0*X.~| 


(32) 


where 


k 

a 

3 


f\  ♦  e“ 


-1 

tan  e 


dw 

9t 


♦  w 


dw 

ds 


(33a) 

(33b) 

(33c) 


The  finite-difference  approximations  to  Eq.  (32)  are  written  along  the  stream¬ 
line  direction  (see  Pig.  2)  at 
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$ 


I 


::: 


'fa 


y/.j 


,yl 

1 

•*T 


V,‘l 


-  h-V  ♦ 


♦  i  <gi-i/2  * 


(q",n~l  -  ^:rl>  ♦  *p 


m.n-l..p 

9J-l/2)ej-l/2 


1  (jLi.n  v«.n-l 

2  ^  3-1/2  3-1/2 


m.n-1 


4,,)-l/2 


where 


*3- 1/2 


Vco“m/2 


The  relation  between  8 
end  ( 1-3/2, n-1/2)  ere 


j-1/2 


end  those  values  of  8  centered  at  ( 1- 1/2, n- 1/2) 


.1-3/2  .1-1/2 

93-l/2  ~  .1-1/2 


zll2 — (sp  -  S  )  ♦  el_3/2 
*1-3/2  *1-1/2  1-3/2  1-1/2 


3-1/2 


The  solution  of  Sq.  (12)  subject  to  the  boundary  conditions  given  by  Bq.  (13) 

Is  achieved  by  solving  the  system  of  equations  given  by  Bqs.  (18),  (19)  and 
(20)  (standard  scheme)  when  there  Is  no  flow  reversal.  When  calculated  results 
reveal  flow  reversal  (e^  <  0),  further  iterations  at  that  location  make  use 
of  the  characteristic  scheme  which  seeks  the  solution  of  Bqs.  (18),  (19a,b,c) 


and  (32)  for  e^  <  0  and  the  regular  seb 


for  e  >  0.  This  switch  from 


one  scheme  to  another  continues  to  allow  quadratic  convergence  and  ensures  that 
numerical  Instabilities  are  avoided  provided  that  the  step  lengths  in  the  t- 
and  s-dlrectlons  are  "properly"  selected  as  we  shall  discuss  in  the  following 
section. 
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4.0  THE  QUESTION  OP  SINGULARITY  ON  AN  OSCILLATING  AIRPOIL 


The  problem  of  a  circular  cylinder  impulsively  started  from  rest  has  served  as 
a  model  problem  with  which  to  examine  unsteady  boundary  layers  and  the  nature 
of  their  solutions  in  the  presence  of  large  flow  reversal.  Noteworthy  contri¬ 
butions  have  been  made  by  Cebeci  (1979),  van  Dommelen  and  Shen  (1982),  Cowley 
(1983)  and  Ingham  (1984),  and  show  that  at  large  times  the  distribution  of 
displacement  thickness  has  a  steep  rise  near  the  location  of  zero  wall  shear 
and  with  consequent  tendency  for  calculations  to  break  down.  The  suggested 
values  for  the  time  and  location  of  zero  wall  shear  and  peaking  of  the  dis¬ 
placement  velocity,  [d/ds(u  4*)]  vary  slightly,  probably  due  to  differences 
in  the  calculation  methods.  The  explanation  for  the  breakdown  of  the  calcula¬ 
tions  has  been  provided  by  Cebeci  (1986)  who  demonstrated  that  numerical  calc¬ 
ulations  must  satisfy  a  CPL-llke  stability  criterion,  if  this  is  done,  it  is 
expected  that  the  location  of  the  singularity  associated  with  unsteady  flow  and 
large  time  will  correspond  exactly  to  that  of  steady  flow,  namely  6  »  105°, 
rather  than  0  *  111*.  The  same  situation  cannot  be  expected  with  oscillating 
airfoils  where  the  solutions  are  cyclic  and  do  not  tend  to  a  steady  state. 

The  present  study  examines  the  nature  of  solutions  to  the  boundary- layer  equa¬ 
tions  for  the  flow  on  an  oscillating  airfoil,  which  can  give  rise  to  extensive 
regions  of  flow  reversal  and  separation.  Here  flow  reversal  implies  the  exist¬ 
ence  of  negative  wall  shear  and  separation  is  taken  to  correspond  to  situations 
where  calculations  with  a  prescribed  pressure  distribution  breakdown  due  to  a 

singularity.  The  calculations  were  made  for  three  values  of  w  with  5  *  1  and 

o 

A  *  -1/2.  With  the  choice  of  u  >  0.001,  0.01  and  0.10,  the  maximum  value  of 
defined  by 

5  cc  =*  l  (1  +  A  sinwt) 
ett  o 

is  sufficient  to  provoke  separation  with  a  strong  singularity.  Por  example. 

the  maximum  value  of  £  __  is  1.5  at  wt  *  270  and  the  flow  conditions 

eff 

closely  resemble  a  steady  separated  flow  at  the  smaller  frequencies  u>  =  0.001 
and  0.01.  Since  the  value  of  l  ..  corresponding  to  steady  flow  separation 

6ll 

is  1.115,  we  would  expect  the  calculations  to  break  down  before  wt  -  270°  due 
to  the  singularity.  Por  the  higher  frequency  case  (w  =  0.10),  we  expect  the 
solutions  to  break  down  later  than  wt  3  270°  with  flow  reversals  occurring  in 
the  range  270°  <  wt  <  360°. 
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The  calculations  war#  arranged  to  parallel  those  previously  performed  for  a 

circular  cylinder  and  reported  by  Cebecl  (1986).  Thus  both  the  zig-zag  and 

the  characteristic  box  schemes  ware  used  first  with  t lme  and  distance  steps 

which  were  chosen  arbitrarily  and  subsequently  with  values  in  agreement  with 

the  stability  criterion.  The  results  of  Pig.  3  for  w  -  0.10  were  obtained 

with  the  zig-zag  box  scheme  by  Cebecl.  Khattab  and  Schimke  (1984)  for  a  At 

spacing  specified  such  that  A£^  -  0.01  up  to  t  *  1.7,  At^  *  0.005  for  1.7 

<  t  <  4  and  At  *  0.01  for  4  <  t  <  8;  the  time  steps  k  were  10  degrees 
1  n 

degrees  for  0  <  wt  <  260.  5  degrees  for  260  <  wt  <  295,  and  1.25  degrees 

for  295  <  wt  <  360.  The  calculations  broke  down  at  wt  ■  310*.  indicating 

flow  separation  at  this  location. 

Pig.  3a  shows  that  the  variation  of  the  displacement  thickness 


4* 


1_L_L 


L. 

Cf 


(25) 


is  generally  saooth  except  in  the  neighborhood  of  t  *  2.12  and  for  wt  * 

308.75*.  The  first  sign  of  irregularity  is  the  steepening  of  the  slope  of 

4*  when  wt  *  300*  and  a  local  maximum  of  4*  occurs  at  i  »  2.12  when 

wt  *  308.75*.  When  the  sane  results  are  plotted  for  a  displacenent  velocity, 

(d/dU(u  4*),  (Pig.  3b),  we  observe  that  the  steepening  of  the  displacenent 
© 

velocity  as  the  peak  moves  form  \  *  2.125  to  2.08  with  wt  changing  fron  300 
to  308.75  degrees.  It  should  be  noted  that  the  maximal  value  of  displaceewnt 
velocity  moves  towards  the  separation  point  with  increasing  wt  and  the  sane 
behavior  will  be  shown  to  occur  for  the  circular  cylinder  discussed  below,  as 

m 

shown  in  Pig.  3c,  the  wall  shear  parameter  f^  shows  no  signs  of  irregu¬ 
larity  for  wt  <  308.75*  but  a  deep  mlnimm  in  f*  occurs  near  l  -  2.15, 
l.e.  near  the  peak  of  4*. 


It  is  Interesting  and  useful  to  compare  the  results  presented  in  Pig.  3  for  an 
oscillating  airfoil  with  those  obtained  by  Cebecl  (1982)  for  a  circular  cylin¬ 
der  started  impulsively  from  rest.  Pig.  4,  and  obtained  with  the  sane  zig-zag 
scheme.  As  in  the  case  of  the  oscillating  airfoil,  the  flow  separates  and 
remains  smooth  up  to  the  separation  point.  However,  just  downstream  of  sepa¬ 
ration  with  increasing  time,  a  singularity  seems  to  develop  in  the  neighborhood 
of  0  *  112*  and  t  *  3.0  and  it  was  not  possible  to  continue  the  boundary- 
layer  calculations  beyond  this  time  and  angular  location.  Prom  Pig.  4a  we  see 
that  the  variation  of  displacement  thickness  is  smooth  for  values  of  6  less 
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than  108*  and  it  bag ins  to  steepen  thereafter.  The  sane  results  are  plotted 
in  Fig.  4b  to  danonstrate  that  the  displacement  velocity  exhibits  a  maximum 
which  Increases  rapidly  with  time,  as  in  Pig.  3b,  with  the  maximum  shifting 
towards  the  location  of  separation  with  increasing  time.  The  results  of  locai 
skin-friction  coefficients,  Pig.  4c,  follow  similar  trends  to  those  obtained 
for  the  oscillating  airfoil  with  the  distributions  passing  through  zero  with 
no  signs  of  Irregularity  and  no  breakdown  before  the  time  corresponding  to  the 
singularity. 

The  calculations  which  led  to  Pig.  3  were  repeated  with  the  characteristic  box 

scheaM  using  the  sane  coarse  variations  of  k  and  A£.  and  the  results 

n  1 

were  identical  to  those  obtained  with  the  zig-zag  scheme  up  to  wt  *  280.  At 
wt  ■  282.5,  the  solutions  of  the  zig-zag  scheme  were  smooth  and  free  of  wig¬ 
gles  but  those  of  the  characteristic  box  schesw  exhibited  oscillations  in 
f  which  led  to  their  breakdown.  The  solutions  with  the  zig-zag  scheme, 
however,  continued  without  numerical  difficulties  until  wt  *  310,  where 
oscillations  appeared  and  led  to  the  breakdown  of  the  solutions  at  the  next 
tine  step. 

The  characteristic  box  was  used  subsequently  with  values  of  Afc^  fixed  as 

before  and  with  values  of  k  determined  in  accord  with  the  stability  requlre- 

n 

sent  as  shown  in  Table  1.  This  procedure  avoided  the  breakdown  of  the  solu¬ 
tions  and,  as  can  be  seen  from  Pig.  5,  the  maximum  value  of  8  Increases  con¬ 
siderably  with  wt  so  that  the  solutions  required  correspondingly  smaller 
values  of  the  time  step.  It  is  Interesting  to  note  that  the  wall-shear  distri¬ 
butions  of  Pig.  6  are  uninfluenced  by  the  mesh  at  wt  *  280  and  310  but,  for 
wt  >  310.  the  coarse  mesh  leads  to  large  values  of  8  and  breakdown  of  the 
solutions. 

Plgure  7a  shows  the  distributions  of  displacement  thickness  for  values  of  wt 
from  260*  up  to  360*  and  completes  the  cycle.  The  results  up  to  300*  were 
identical  with  those  of  Pig.  3a  with  rapid  increase  of  the  displacement  thick¬ 
ness  corresponding  to  the  increasing  extent  of  flow  reversal,  as  shown  by  the 
wall-shear  distributions  of  Pig.  7b.  It  can  also  be  seen  from  this  figure  that 
the  maximum  displacement  thickness  and  minimum  wall  shear  move  upstream  with 
Increasing  wt  for  values  of  wt  up  to  324.5*;  this  feature  was  also  observed 
in  the  calculations  performed  for  the  circular  cylinder  and  shown  in  Pig.  4. 
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Table  1.  The  distribution  of  step  sizes  in  wt  for  wt  »  0.1  in  accordance 
with  the  requirements  of  the  stability  parameter  0 

k 

wt  n 


0  - 

240 

O 

O 

240  - 

255 

5* 

255  - 

261 

3° 

261  - 

265 

2° 

265  - 

284 

1° 

284  - 

305 

0.5° 

305  - 

320 

0.25 

320  - 

360 

0.5* 

The  results  obtained  with  the  zig-zag  scheme  and  values  of  k^  determined  by 
the  characteristic  schesm  for  the  oscillating  airfoil  were  identical  to  those 
discussed  above,  and  siaillar  correspondence  was  obtained  with  the  calculations 
performed  for  the  circular  cylinder. 

Figures  8  and  9  show  the  distributions  of  wall  shear  and  displacement  thickness 
for  two  ssuiller  frequencies  w  -  0.01  and  0.001.  As  expected,  the  critical 
value  of  the  reduced  angle  which  corresponds  to  separation,  is  smaller  than 
that  for  the  higher  frequency  and  closer  to  that  of  the  steady  state,  l  * 

1.16.  For  w  *  0.01,  the  breakdown  of  the  solutions  occurs  at  wt  -  226*. 
which  corresponds  to  an  effective  reduced  angle  of  *  1.360;  for  w  -  0.001, 

the  corresponding  values  are  wt  *  204*  and  1.203.  Ve  also  note  from  Figures  8a 
and  8b  that  the  flow  is  a  "little”  unsteady  even  at  these  frequencies,  and  the 
solutions  do  not  break  down  with  the  appearance  of  flow  reversal,  which 
Increases  in  extent  as  w  changes  from  0.001  and  0.01. 
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5.0  INTERACTION  AS  AN  ANSWER  TO  THE  QUESTION  OF  SINGULARITY 

The  interaction  procedure  discussed  in  Section  3  has  been  applied  to  the  flow 
problem  examined  previously  in  Section  4  with  the  standard  method  and  the 
results  are  shown  in  Figures  10  to  14  and  discussed  below.  In  contrast  to  the 
standard  problem,  which  makes  the  implicit  assumption  of  infinite  Reynolds 
number,  the  interaction  requires  specification  of  a  finite  Reynolds  number.  A 
thickness  ratio  t  has  also  to  be  specified  and,  since  the  definition  of  e 
involves  R  and  t,  the  calculations  are  performed  for  specified  values  of  e. 

In  all  cases  shown,  the  calculations  made  use  of  time  steps  determined  by  the 
characteristic  scheme  in  agreement  with  the  stability  requirement.  This  was 
not  done  in  the  calculations  of  Cebeci  et  al.  (1984)  and  the  solutions  exhib¬ 
ited  oscillations  which  steamed  from  the  numerical  method. 

The  present  calculations  were  performed  in  the  following  way.  For  all  values 
of  time  with  ut  ranging  from  0  to  360*,  the  standard  method  and  the  leading- 
edge  region  procedure  of  Cebeci  and  Carr  (1981)  were  used  to  generate  initial 
conditions  at  a  short  distance  from  the  leading  edge,  £  *  0.5.  with  these 
Initial  conditions  and  for  each  value  of  ut,  the  inverse  method  was  used  to 
calculate  the  unsteady  flow  from  £  -  0.5  to  10,  for  the  specified  value  of 
c.  since  the  system  of  equations  is  now  elliptic,  sweeps  in  the  ^direc¬ 
tion  were  necessary  to  achieve  a  converged  solution;  around  three  sweeps  were 
required  where  flow  reversal  was  encountered  and  a  single  sweep  sufficed  where 
it  did  not.  It  is  to  be  expected  that  the  value  of  c  will  influence  the 
number  of  sweeps  and,  since  it  is  linked  to  physical  parameters,  will  affect 
the  singularity  and  the  size  of  the  bubble. 

4 

Figures  10  and  11  show  the  results  for  w  -  0.001  and  0.01  with  c  ■  10  .  They 
are  nearly  the  same  as  those  obtained  by  the  standard  method  and  shown  in  Fig¬ 
ures  8  and  9  prior  to  flow  reversal  where  the  influence  of  the  Reynolds  number 
is  small  and  increase  after  flow  reversal.  In  the  case  of  u  3  0.001,  for 
example,  the  standard  method  predicts  flow  reversal  around  £  =  1.19  (see 

6t  L 

Pig.  9)  and  with  interaction  (Fig.  10)  this  effective  angle  is  between  1.219 

m 

and  1.254.  The  maximum  negative  value  of  the  wall-shear  parameter  f^  obtained 

with  the  standard  method  is  around  -0.03  at  £  -  1.199  and  may  be  compared 

„  eff 

with  the  maximum  value  of  f^  of  -0.14  at  £fiff  =  1.286  obtained  with  inter¬ 
action.  As  expected,  the  interaction  allows  the  calculations  to  be  performed 
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at  higher  angles  of  attack  than  those  achieved  with  the  standard  method.  For 
oj  =  0.001,  the  maximum  a  for  which  calculations  can  be  performed  with 

6ll 

the  standard  method  is  1.199  with  breakdown  occurring  at  l  __  =  1.209;  the 

ef  f 

corresponding  values  with  interaction  are  1.286  and  1.287.  Comparison  of  wall- 

shear  results  with  both  procedures  and  u  =  0.001  indicates  that  the  extent 

of  the  recirculation  region  is  around  0.5  for  the  standard  case,  and 

around  2.5  for  the  interactive  case.  The  solutions  do  not  have  a  singularity 

in  the  former  case  but  do  contain  flow  reversals  and  this  suggests  that  time- 

dependent  flows  can  be  calculated  without  using  an  Inverse  procedure.  As  the 

angle  of  attack  exceeds  %  ..  =  1.199  for  u  =  0.001,  a  singularity  devel- 

ett 

ops  and  requires  an  inverse  procedure  as  in  two-dimensional  steady  flows.  This 
procedure  allows  the  calculation  of  larger  regions  of  reverse  flow  where  the 
flow  is  now  separated. 

We  see  a  similar  picture  with  the  greater  unsteadiness  corresponding  to  u  = 

0.01,  for  which  the  standard  method  allows  calculations  up  to  an  effective 

angle  of  attack  of  1.354  (Fig.  8a),  a  value  considerably  higher  than  1.199 

obtained  at  w  *  0.001.  The  first  flow  reversal  occurs  shortly  after  = 

1.294  and  breakdown  occurs  at  5  „  *  1.360  with  maximum  negative  wall  shear 

eff 

values  of  -0.14  at  E  „  *  1.354  and  -0.035  at  5  *  1.315.  The  extent  of 

eff  eff 

the  maximum  reverse-flow  region  is  now  1.5,  considerably  larger  than  for  u>  = 
0.001,  and  indicates  that  the  more  unsteady  nature  of  the  flow  produces  a  big¬ 
ger  region  of  reverse  flow  free  from  singularities.  For  this  value  of  w,  the 
interactive  scheme  Increases  the  value  of  l  __  for  which  solutions  can  be 
obtained  to  1.424  with  breakdown  occurring  shortly  after  this  value  at  1.428 
(see  Fig.  11).  The  first  flow  reversal  occurs  after  £  *  1.315  with 

maximum  negative  wall  shear  equal  to  -0.19  at  l  __  *  1.424,  and  the  extent  of 

eff 

the  recirculation  region  has  now  increased  by  about  30%.  Comparison  of  maxi- 

n 

mum  wall  shear  values,  fw,  at  a  similar  value  of  indicates  that 

those  computed  with  the  interactive  scheme  are  lower  than  those  with  the  stand- 

IV 

ard  scheme  so  that,  for  example,  the  interactive  scheme  gives  (f  )  -  -  0.04 

’  w  max 

at  E  =  1.36  compared  to  -0.14  at- 5  __  =  1.354  with  the  standard  method 
eff  eff 

(Pig.  8a). 

Figure  12  shows  that  the  size  of  the  reverse- flow  region  Increases  with 
Reynolds  number  but  the  effective  angle  of  attack  for  which  solutions  can  be 

4 

obtained  is  only  slightly  reduced,  changing  from  1.428  for  c  =*  10  to 
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5 

around  1.415  for  e  =  10  .  It  is  interesting  to  note  that  the  interactive 

4 

solutions  do  not  have  any  flow  reversal  at  =  1.315  with  c  =  10  . 

4 

Figures  13  and  14  show  the  results  for  u  =  0.1  with  values  of  c  of  10 
and  10  and  they  are  again  similar  to  those  obtained  by  the  standard  method, 
as  shown  in  Figure  7,  prior  to  flow  reversal  where  the  influence  of  Reynolds 
number  is  small.  After  flow  reversal,  the  differences  between  the  results 
obtained  with  the  standard  and  interactive  methods  increase  as  the  Reynolds 
number  decreases.  It  is  clear  that  the  solutions  are  free  from  the  numerical 
"wiggles"  encountered  when  the  stability  criterion  was  not  met. 

Comparison  of  results  obtained  at  the  two  Reynolds  numbers  for  w  =  0.1  indi¬ 
cates  that  the  interaction  does  not  reduce  the  maximum  negative  value  of  the 

it 

wall  shear  parameter  as  it  did  with  lower  frequencies.  For  example,  f 

at  ut  =  360°  is  around  -0.19  with  the  standard  scheme  and  around  -0.30  at 

4  5 

e  *  10  and  around  -0.35  at  c  =  10  with  the  interactive  method.  The 
maximum  value  of  negative  wall  shear  calculated  with  interaction  is  consider¬ 
ably  greater  them  its  corresponding  value  obtained  with  the  standard  method  at 
the  end  of  one  complete  cycle.  Furthermore,  the  behavior  of  the  wall  shear  is 

N 

not  monotonic  without  interaction  so  that,  for  example,  fw  reaches  a  max¬ 
imum  value  equal  to  -0.25  around  wt  *  331°  and  then  decreases  to  -0.195  at 
wt  =  360°.  With  interaction  this  is  not  the  case  with  the  maximum  negative 

m 

value  of  f  continuously  Increasing  with  u>t. 

The  results  of  Figures  7,  13  and  14  for  u  *  0.1  are  for  an  unsteady  flow  and 
are  unlike  those  for  two  other  values  of  w  in  that  they  are  free  from  singu¬ 
larities.  For  this  reason,  even  though  the  results  in  the  reverse  flow  region 
and  thereafter  are  different  due  to  the  Reynolds  number  effect,  the  extent  of 
the  reverse-flow  region  is  essentially  the  same  and  is  consistent  with  the 
results  obtained  at  lower  frequencies  in  the  absence  of  flow  separation  even 
though  the  extent  of  the  reverse- flow  region  is  reduced  at  the  lower  Reynolds 
numbers. 

The  results  obtained  with  u  =  0.001  can  usefully  be  compared  with  the  steady- 
state  results  of  Cebecl  et  al.  (1981)  shown  in  Figure  15.  We  might  expect  that 
the  small  unsteadiness  associated  with  this  frequency  will  lead  to  results  very 
similar  to  those  of  steady  state.  Inspection  of  Figures  10  and  15  shows  that 
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although  this  is  correct  in  general  terras,  the  answers  are  quantitatively  dif¬ 
ferent.  As  can  be  seen,  the  maximum  effective  angle  at  which  solutions  can  be 
obtained  is  greater  in  the  unsteady  case  by  some  7%.  There  are  differences  in 
the  two  calculation  procedures  but  it  is  unlikely  that  they  are  responsible  for 
this  difference.  On  the  other  hand,  it  is  possible  that  the  difference  in  the 
negative  wall  shear  values  may  have  been  influenced  by  the  use  of  the  FLARE 
approximation  in  the  steady-state  solutions.  Nevertheless,  the  unsteady  nature 
of  the  flow  with  w  =  0.001  is  clear,  in  spite  of  this  very  low  reduced 
f  requency . 
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6 . 0  CONCLUDING  REMARKS 


The  following  principle  conclusions  may  be  drawn  from  the  preceding  text. 

1.  A  calculation  method  has  been  developed  to  represent  flows  around  oscil¬ 
lating  airfoils.  It  is  based  on  a  similar  approach  used  for  steady  flows 
with  separation  and  involves  interaction  between  inviscid  and  viscous-flow 
equations.  The  coupling  technique  is  similar  to  that  described  by  Veldman 
(1981)  and  Cebeci  et  al.  (1981)  for  steady  flows.  This  interactive  method 
has  been  used  to  calculate  separation  and  reattachment  near  the  leading 
edge  of  a  thin  oscillating  airfoil  and  has  been  shown  to  give  rise  to  rapid 
convergence  similar  to  that  obtained  in  steady  flows,  Cebeci  et  al.  (1986). 

2.  The  accuracy  of  the  results  obtained  from  the  solution  of  the  boundary- 
layer  equations  has  been  examined  with  emphasis  on  regions  of  flow  reversal 
and  separation  where  the  characteristic  box  scheme  is  used.  Attempts  to 
improve  accuracy  by  ad  hoc  changes  to  the  finite-difference  mesh  failed  and 
revealed  the  need  for  a  procedure  which  would  automatically  guarantee 
accuracy  by  the  selection  of  an  appropriate  mesh.  This  was  achieved 
through  a  stability  criterion,  similar  to  that  of  Courant,  Friedrichs  and 
Lewy.  The  combination  of  this  requirement  and  the  characteristic  box 
scheme  led  to  accurate  solutions  and  showed  that  the  mesh  requirements 
were  extremely  severe  in  the  region  of  large  flow  reversals. 

3.  Calculations  have  been  performed  for  a  range  of  reduced  frequencies  from  0 
to  0.1.  They  show  that  increased  unsteadiness  allows  results  to  be 
obtained  at  higher  angles  of  attack  before  the  solutions  break  down;  indeed 
in  the  case  of  the  highest  frequency  there  was  no  breakdown.  The  calcula¬ 
tions  with  the  standard  method  led  to  regions  of  flow  reversal  which  were 
limited  in  their  extent  by  the  singularity  except  at  the  highest  frequency. 
The  interactive  procedure  removed  this  singularity  and  resulted  in  larger 
regions  of  flow  reversal  which  involved  separation  at  higher  angles  of 
attack. 

The  calculated  maximum  angles  of  attack  were,  however,  modest  and  regions 
of  separated  flow  were  small.  This  is  consistent  with  the  behavior  of 
steady  laminar  flows  which  can  only  sustain  small  separation  bubbles. 
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The  unsteady  nature  of  the  flow  at  the  highest  frequency  allowed  the  calc¬ 
ulation  of  large  regions  of  flow  reversal  and  it  is  expected  that  yet 
higher  frequencies  will  lead  to  even  larger  regions  of  flow  reversal.  This 
in  turn  will  permit  calculations  to  be  performed  at  larger  angles  of  attack 
where  the  occurrence  of  the  singularity  will  require  the  use  of  the  inter¬ 
active  procedure.  The  gains  in  angles  of  attack  are  again  likely  to  be 
limited  by  the  ability  of  the  laminar  flow  to  sustain  separation  bubbles. 


The  interactive  scheme,  incorporating  the  solution  of  the  boundary- layer 
equations  by  the  characteristic  box  scheme  and  with  the  numerical  mesh 
determined  in  accordance  with  the  stability  criterion,  has  been  used  to 
calculate  the  laminar  flow  for  a  model  problem.  The  numerical  aspects  of 
this  procedure  have  been  thoroughly  tested  and  shown  to  be  general  so  that 
it  can  be  used  for  the  solution  of  laminar  and  turbulent  flows  over  air¬ 
foils  of  practical  relevance. 
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Effect  of  interaction  on  the  variation  of  (a)  wall  shear  f"  and  (b)  displacement  thickness 
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Figure  13.  Effect 


